{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<style>\n",
       "@import url('http://fonts.googleapis.com/css?family=Source+Code+Pro');\n",
       "@import url('http://fonts.googleapis.com/css?family=Vollkorn');\n",
       "@import url('http://fonts.googleapis.com/css?family=Arimo');\n",
       "@import url('http://fonts.googleapis.com/css?family=Fira_sans');\n",
       "\n",
       "    div.cell{\n",
       "        width: 900px;\n",
       "        margin-left: 0% !important;\n",
       "        margin-right: auto;\n",
       "    }\n",
       "    div.text_cell code {\n",
       "        background: transparent;\n",
       "        color: #000000;\n",
       "        font-weight: 600;\n",
       "        font-size: 11pt;\n",
       "        font-style: bold;\n",
       "        font-family:  'Source Code Pro', Consolas, monocco, monospace;\n",
       "   }\n",
       "    h1 {\n",
       "        font-family: 'Open sans',verdana,arial,sans-serif;\n",
       "\t}\n",
       "\t\n",
       "    div.input_area {\n",
       "        background: #F6F6F9;\n",
       "        border: 1px solid #586e75;\n",
       "    }\n",
       "\n",
       "    .text_cell_render h1 {\n",
       "        font-weight: 200;\n",
       "        font-size: 30pt;\n",
       "        line-height: 100%;\n",
       "        color:#c76c0c;\n",
       "        margin-bottom: 0.5em;\n",
       "        margin-top: 1em;\n",
       "        display: block;\n",
       "        white-space: wrap;\n",
       "        text-align: left;\n",
       "    } \n",
       "    h2 {\n",
       "        font-family: 'Open sans',verdana,arial,sans-serif;\n",
       "        text-align: left;\n",
       "    }\n",
       "    .text_cell_render h2 {\n",
       "        font-weight: 200;\n",
       "        font-size: 16pt;\n",
       "        font-style: italic;\n",
       "        line-height: 100%;\n",
       "        color:#c76c0c;\n",
       "        margin-bottom: 0.5em;\n",
       "        margin-top: 1.5em;\n",
       "        display: block;\n",
       "        white-space: wrap;\n",
       "        text-align: left;\n",
       "    } \n",
       "    h3 {\n",
       "        font-family: 'Open sans',verdana,arial,sans-serif;\n",
       "    }\n",
       "    .text_cell_render h3 {\n",
       "        font-weight: 200;\n",
       "        font-size: 14pt;\n",
       "        line-height: 100%;\n",
       "        color:#d77c0c;\n",
       "        margin-bottom: 0.5em;\n",
       "        margin-top: 2em;\n",
       "        display: block;\n",
       "        white-space: wrap;\n",
       "        text-align: left;\n",
       "    }\n",
       "    h4 {\n",
       "        font-family: 'Open sans',verdana,arial,sans-serif;\n",
       "    }\n",
       "    .text_cell_render h4 {\n",
       "        font-weight: 100;\n",
       "        font-size: 14pt;\n",
       "        color:#d77c0c;\n",
       "        margin-bottom: 0.5em;\n",
       "        margin-top: 0.5em;\n",
       "        display: block;\n",
       "        white-space: nowrap;\n",
       "    }\n",
       "    h5 {\n",
       "        font-family: 'Open sans',verdana,arial,sans-serif;\n",
       "    }\n",
       "    .text_cell_render h5 {\n",
       "        font-weight: 200;\n",
       "        font-style: normal;\n",
       "        color: #1d3b84;\n",
       "        font-size: 16pt;\n",
       "        margin-bottom: 0em;\n",
       "        margin-top: 0.5em;\n",
       "        display: block;\n",
       "        white-space: nowrap;\n",
       "    }\n",
       "    div.text_cell_render{\n",
       "        font-family: 'Fira sans', verdana,arial,sans-serif;\n",
       "        line-height: 125%;\n",
       "        font-size: 115%;\n",
       "        text-align:justify;\n",
       "        text-justify:inter-word;\n",
       "    }\n",
       "    div.output_subarea.output_text.output_pyout {\n",
       "        overflow-x: auto;\n",
       "        overflow-y: scroll;\n",
       "        max-height: 50000px;\n",
       "    }\n",
       "    div.output_subarea.output_stream.output_stdout.output_text {\n",
       "        overflow-x: auto;\n",
       "        overflow-y: scroll;\n",
       "        max-height: 50000px;\n",
       "    }\n",
       "    div.output_wrapper{\n",
       "        margin-top:0.2em;\n",
       "        margin-bottom:0.2em;\n",
       "}\n",
       "\n",
       "    code{\n",
       "      font-size: 70%;\n",
       "    }\n",
       "    .rendered_html code{\n",
       "    background-color: transparent;\n",
       "    }\n",
       "    ul{\n",
       "        margin: 2em;\n",
       "    }\n",
       "    ul li{\n",
       "        padding-left: 0.5em; \n",
       "        margin-bottom: 0.5em; \n",
       "        margin-top: 0.5em; \n",
       "    }\n",
       "    ul li li{\n",
       "        padding-left: 0.2em; \n",
       "        margin-bottom: 0.2em; \n",
       "        margin-top: 0.2em; \n",
       "    }\n",
       "    ol{\n",
       "        margin: 2em;\n",
       "    }\n",
       "    ol li{\n",
       "        padding-left: 0.5em; \n",
       "        margin-bottom: 0.5em; \n",
       "        margin-top: 0.5em; \n",
       "    }\n",
       "    ul li{\n",
       "        padding-left: 0.5em; \n",
       "        margin-bottom: 0.5em; \n",
       "        margin-top: 0.2em; \n",
       "    }\n",
       "    a:link{\n",
       "       font-weight: bold;\n",
       "       color:#447adb;\n",
       "    }\n",
       "    a:visited{\n",
       "       font-weight: bold;\n",
       "       color: #1d3b84;\n",
       "    }\n",
       "    a:hover{\n",
       "       font-weight: bold;\n",
       "       color: #1d3b84;\n",
       "    }\n",
       "    a:focus{\n",
       "       font-weight: bold;\n",
       "       color:#447adb;\n",
       "    }\n",
       "    a:active{\n",
       "       font-weight: bold;\n",
       "       color:#447adb;\n",
       "    }\n",
       "    .rendered_html :link {\n",
       "       text-decoration: underline; \n",
       "    }\n",
       "    .rendered_html :hover {\n",
       "       text-decoration: none; \n",
       "    }\n",
       "    .rendered_html :visited {\n",
       "      text-decoration: none;\n",
       "    }\n",
       "    .rendered_html :focus {\n",
       "      text-decoration: none;\n",
       "    }\n",
       "    .rendered_html :active {\n",
       "      text-decoration: none;\n",
       "    }\n",
       "    .warning{\n",
       "        color: rgb( 240, 20, 20 )\n",
       "    } \n",
       "    hr {\n",
       "      color: #f3f3f3;\n",
       "      background-color: #f3f3f3;\n",
       "      height: 1px;\n",
       "    }\n",
       "    blockquote{\n",
       "      display:block;\n",
       "      background: #fcfcfc;\n",
       "      border-left: 5px solid #c76c0c;\n",
       "      font-family: 'Open sans',verdana,arial,sans-serif;\n",
       "      width:680px;\n",
       "      padding: 10px 10px 10px 10px;\n",
       "      text-align:justify;\n",
       "      text-justify:inter-word;\n",
       "      }\n",
       "      blockquote p {\n",
       "        margin-bottom: 0;\n",
       "        line-height: 125%;\n",
       "        font-size: 100%;\n",
       "      }\n",
       "</style>\n",
       "<script>\n",
       "    MathJax.Hub.Config({\n",
       "                        TeX: {\n",
       "                           extensions: [\"AMSmath.js\"]\n",
       "                           },\n",
       "                tex2jax: {\n",
       "                    inlineMath: [ ['$','$'], [\"\\\\(\",\"\\\\)\"] ],\n",
       "                    displayMath: [ ['$$','$$'], [\"\\\\[\",\"\\\\]\"] ]\n",
       "                },\n",
       "                displayAlign: 'center', // Change this to 'center' to center equations.\n",
       "                \"HTML-CSS\": {\n",
       "                    scale:100,\n",
       "                        availableFonts: [],\n",
       "                        preferredFont:null,\n",
       "                        webFont: \"TeX\",\n",
       "                    styles: {'.MathJax_Display': {\"margin\": 4}}\n",
       "                }\n",
       "        });\n",
       "</script>\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# format the book\n",
    "%matplotlib inline\n",
    "import sys\n",
    "from __future__ import division, print_function\n",
    "import sys\n",
    "sys.path.insert(0,'../code')\n",
    "import book_format\n",
    "book_format.load_style('../code')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Two Dimensions\n",
    "\n",
    "\n",
    "## Paintball\n",
    "\n",
    "Paintball is a sport in which competing teams try to shoot each other\n",
    "with guns that fire paint-filled pellets that break on impact, leaving a\n",
    "colorful mark on the target. It is usually played in an arena decorated\n",
    "with barriers and other objects that can be used as cover.\n",
    "\n",
    "Suppose you are playing paintball in an indoor arena 30 feet wide and 50\n",
    "feet long. You are standing near one of the 30 foot walls, and you\n",
    "suspect that one of your opponents has taken cover nearby. Along the\n",
    "wall, you see several paint spatters, all the same color, that you think\n",
    "your opponent fired recently.\n",
    "\n",
    "The spatters are at 15, 16, 18, and 21 feet, measured from the\n",
    "lower-left corner of the room. Based on these data, where do you think\n",
    "your opponent is hiding?\n",
    "\n",
    "Figure [fig.paintball] shows a diagram of the arena. Using the\n",
    "lower-left corner of the room as the origin, I denote the unknown\n",
    "location of the shooter with coordinates $\\alpha$ and $\\beta$, or\n",
    "`alpha` and `beta`. The location of a spatter is\n",
    "labeled `x`. The angle the opponent shoots at is $\\theta$ or\n",
    "`theta`.\n",
    "\n",
    "The Paintball problem is a modified version of the Lighthouse problem, a\n",
    "common example of Bayesian analysis. My notation follows the\n",
    "presentation of the problem in D.S. Sivia’s, *Data Analysis: a\n",
    "Bayesian Tutorial, Second Edition* (Oxford, 2006).\n",
    "\n",
    "You can download the code in this chapter from\n",
    "<http://thinkbayes.com/paintball.py>. For more information see\n",
    "Section [download]."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The suite\n",
    "\n",
    "![Diagram of the layout for the paintball problem.](figs/paintball.pdf)\n",
    "\n",
    "[fig.paintball]\n",
    "\n",
    "To get started, we need a Suite that represents a set of hypotheses\n",
    "about the location of the opponent. Each hypothesis is a pair of\n",
    "coordinates: `(alpha, beta)`.\n",
    "\n",
    "Here is the definition of the Paintball suite:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import thinkbayes\n",
    "\n",
    "class Paintball(thinkbayes.Suite, thinkbayes.Joint):\n",
    "\n",
    "    def __init__(self, alphas, betas, locations):\n",
    "        self.locations = locations\n",
    "        pairs = [(alpha, beta) \n",
    "                 for alpha in alphas \n",
    "                 for beta in betas]\n",
    "        thinkbayes.Suite.__init__(self, pairs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`Paintball` inherits from `Suite`, which we have\n",
    "seen before, and `Joint`, which I will explain soon.\n",
    "\n",
    "`alphas` is the list of possible values for\n",
    "`alpha`; `betas` is the list of values for\n",
    "`beta`. `pairs` is a list of all `(alpha,\n",
    "beta)` pairs.\n",
    "\n",
    "`locations` is a list of possible locations along the wall;\n",
    "it is stored for use in `Likelihood`.\n",
    "\n",
    "![Posterior CDFs for `alpha` and `beta`, given the\n",
    "data.](figs/paintball2.pdf)\n",
    "\n",
    "[fig.paintball2]\n",
    "\n",
    "The room is 30 feet wide and 50 feet long, so here’s the code that\n",
    "creates the suite:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "alphas = range(0, 31)\n",
    "betas = range(1, 51)\n",
    "locations = range(0, 31)\n",
    "\n",
    "suite = Paintball(alphas, betas, locations)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This prior distribution assumes that all locations in the room are\n",
    "equally likely. Given a map of the room, we might choose a more detailed\n",
    "prior, but we’ll start simple."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Trigonometry\n",
    "\n",
    "Now we need a likelihood function, which means we have to figure out the\n",
    "likelihood of hitting any spot along the wall, given the location of the\n",
    "opponent.\n",
    "\n",
    "As a simple model, imagine that the opponent is like a rotating turret,\n",
    "equally likely to shoot in any direction. In that case, he is most\n",
    "likely to hit the wall at location `alpha`, and less likely\n",
    "to hit the wall far away from `alpha`.\n",
    "\n",
    "With a little trigonometry, we can compute the probability of hitting\n",
    "any spot along the wall. Imagine that the shooter fires a shot at angle\n",
    "$\\theta$; the pellet would hit the wall at location $x$, where\n",
    "\n",
    "$$x - \\alpha = \\beta \\tan \\theta$$ Solving this equation for $\\theta$\n",
    "yields $$\\theta = tan^{-1} \\left( \\frac{x - \\alpha}{\\beta} \\right)$$ \n",
    "\n",
    "So given a location on the wall, we can find $\\theta$.\n",
    "\n",
    "Taking the derivative of the first equation with respect to $\\theta$\n",
    "yields\n",
    "\n",
    "$$\\frac{dx}{d\\theta} = \\frac{\\beta}{\\cos^2 \\theta}$$ \n",
    "\n",
    "This derivative is what I’ll call the “strafing speed”, which is the speed of\n",
    "the target location along the wall as $\\theta$ increases. The\n",
    "probability of hitting a given point on the wall is inversely related to\n",
    "strafing speed.\n",
    "\n",
    "If we know the coordinates of the shooter and a location along the wall,\n",
    "we can compute strafing speed:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def StrafingSpeed(alpha, beta, x):\n",
    "    theta = math.atan2(x - alpha, beta)\n",
    "    speed = beta / math.cos(theta)**2\n",
    "    return speed"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`alpha` and `beta` are the coordinates of the\n",
    "shooter; `x` is the location of a spatter. The result is the\n",
    "derivative of `x` with respect to `theta`.\n",
    "\n",
    "![PMF of location given `alpha=10`, for several values of\n",
    "`beta`.](figs/paintball1.pdf)\n",
    "\n",
    "[fig.paintball1]\n",
    "\n",
    "Now we can compute a Pmf that represents the probability of hitting any\n",
    "location on the wall. `MakeLocationPmf` takes\n",
    "`alpha` and `beta`, the coordinates of the\n",
    "shooter, and `locations`, a list of possible values of\n",
    "`x`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def MakeLocationPmf(alpha, beta, locations):\n",
    "    pmf = thinkbayes.Pmf()\n",
    "    for x in locations:\n",
    "        prob = 1.0 / StrafingSpeed(alpha, beta, x)\n",
    "        pmf.Set(x, prob)\n",
    "    pmf.Normalize()\n",
    "    return pmf"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`MakeLocationPmf` computes the probability of hitting each\n",
    "location, which is inversely related to strafing speed. The result is a\n",
    "Pmf of locations and their probabilities.\n",
    "\n",
    "Figure [fig.paintball1] shows the Pmf of location with `alpha =\n",
    "10` and a range of values for `beta`. For all values of\n",
    "beta the most likely spatter location is `x = 10`; as\n",
    "`beta` increases, so does the spread of the Pmf."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Likelihood\n",
    "\n",
    "Now all we need is a likelihood function. We can use\n",
    "`MakeLocationPmf` to compute the likelihood of any value of\n",
    "`x`, given the coordinates of the opponent."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def Likelihood(self, data, hypo):\n",
    "    alpha, beta = hypo\n",
    "    x = data\n",
    "    pmf = MakeLocationPmf(alpha, beta, self.locations)\n",
    "    like = pmf.Prob(x)\n",
    "    return like"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Again, `alpha` and `beta` are the hypothetical\n",
    "coordinates of the shooter, and `x` is the location of an\n",
    "observed spatter.\n",
    "\n",
    "`pmf` contains the probability of each location, given the\n",
    "coordinates of the shooter. From this Pmf, we select the probability of\n",
    "the observed location.\n",
    "\n",
    "And we’re done. To update the suite, we can use `UpdateSet`,\n",
    "which is inherited from `Suite`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "ename": "UnimplementedMethodException",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mUnimplementedMethodException\u001b[0m              Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-7-2c681a88071a>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m()\u001b[0m\n\u001b[1;32m----> 1\u001b[1;33m \u001b[0msuite\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mUpdateSet\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m15\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;36m16\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;36m18\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;36m21\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[1;32mC:\\Users\\rlabbe\\Dropbox\\ThinkBayes\\code\\thinkbayes.py\u001b[0m in \u001b[0;36mUpdateSet\u001b[1;34m(self, dataset)\u001b[0m\n\u001b[0;32m   1195\u001b[0m         \u001b[1;32mfor\u001b[0m \u001b[0mdata\u001b[0m \u001b[1;32min\u001b[0m \u001b[0mdataset\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m   1196\u001b[0m             \u001b[1;32mfor\u001b[0m \u001b[0mhypo\u001b[0m \u001b[1;32min\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mValues\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m-> 1197\u001b[1;33m                 \u001b[0mlike\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mLikelihood\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mdata\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mhypo\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m   1198\u001b[0m                 \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mMult\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mhypo\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mlike\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m   1199\u001b[0m         \u001b[1;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mNormalize\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32mC:\\Users\\rlabbe\\Dropbox\\ThinkBayes\\code\\thinkbayes.py\u001b[0m in \u001b[0;36mLikelihood\u001b[1;34m(self, data, hypo)\u001b[0m\n\u001b[0;32m   1218\u001b[0m         \u001b[0mdata\u001b[0m\u001b[1;33m:\u001b[0m \u001b[0msome\u001b[0m \u001b[0mrepresentation\u001b[0m \u001b[0mof\u001b[0m \u001b[0mthe\u001b[0m \u001b[0mdata\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m   1219\u001b[0m         \"\"\"\n\u001b[1;32m-> 1220\u001b[1;33m         \u001b[1;32mraise\u001b[0m \u001b[0mUnimplementedMethodException\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m   1221\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m   1222\u001b[0m     \u001b[1;32mdef\u001b[0m \u001b[0mLogLikelihood\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mdata\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mhypo\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mUnimplementedMethodException\u001b[0m: "
     ]
    }
   ],
   "source": [
    "suite.UpdateSet([15, 16, 18, 21])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The result is a distribution that maps each `(alpha, beta)`\n",
    "pair to a posterior probability."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Joint distributions\n",
    "\n",
    "When each value in a distribution is a tuple of variables, it is called\n",
    "a **joint distribution** because it represents the\n",
    "distributions of the variables together, that is “jointly”. A joint\n",
    "distribution contains the distributions of the variables, as well\n",
    "information about the relationships among them.\n",
    "\n",
    "Given a joint distribution, we can compute the distributions of each\n",
    "variable independently, which are called the **marginal\n",
    "distributions**.\n",
    "\n",
    "`thinkbayes.Joint` provides a method that computes marginal\n",
    "distributions:\n",
    "\n",
    "```python\n",
    "# class Joint:\n",
    "    def Marginal(self, i):\n",
    "        pmf = Pmf()\n",
    "        for vs, prob in self.Items():\n",
    "            pmf.Incr(vs[i], prob)\n",
    "        return pmf\n",
    "```\n",
    "\n",
    "`i` is the index of the variable we want; in this example\n",
    "`i=0` indicates the distribution of `alpha`, and\n",
    "`i=1` indicates the distribution of `beta`.\n",
    "\n",
    "Here’s the code that extracts the marginal distributions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "marginal_alpha = suite.Marginal(0)\n",
    "marginal_beta = suite.Marginal(1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Figure [fig.paintball2] shows the results (converted to CDFs). The\n",
    "median value for `alpha` is 18, near the center of mass of\n",
    "the observed spatters. For `beta`, the most likely values are\n",
    "close to the wall, but beyond 10 feet the distribution is almost\n",
    "uniform, which indicates that the data do not distinguish strongly\n",
    "between these possible locations.\n",
    "\n",
    "Given the posterior marginals, we can compute credible intervals for\n",
    "each coordinate independently:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "print('alpha CI', marginal_alpha.CredibleInterval(50))\n",
    "print('beta CI', marginal_beta.CredibleInterval(50))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The 50% credible intervals are `(14, 21)` for\n",
    "`alpha` and `(5, 31)` for `beta`. So\n",
    "the data provide evidence that the shooter is in the near side of the\n",
    "room. But it is not strong evidence. The 90% credible intervals cover\n",
    "most of the room!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Conditional distributions {#conditional}\n",
    "\n",
    "![Posterior distributions for `alpha` conditioned on several\n",
    "values of `beta`.](figs/paintball3.pdf)\n",
    "\n",
    "[fig.paintball3]\n",
    "\n",
    "The marginal distributions contain information about the variables\n",
    "independently, but they do not capture the dependence between variables,\n",
    "if any.\n",
    "\n",
    "One way to visualize dependence is by computing **conditional\n",
    "distributions**. `thinkbayes.Joint` provides a method\n",
    "that does that:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def Conditional(self, i, j, val):\n",
    "    pmf = Pmf()\n",
    "    for vs, prob in self.Items():\n",
    "        if vs[j] != val: continue\n",
    "        pmf.Incr(vs[i], prob)\n",
    "\n",
    "    pmf.Normalize()\n",
    "    return pmf"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Again, `i` is the index of the variable we want;\n",
    "`j` is the index of the conditioning variable, and\n",
    "`val` is the conditional value.\n",
    "\n",
    "The result is the distribution of the $i$th variable under the condition\n",
    "that the $j$th variable is `val`.\n",
    "\n",
    "For example, the following code computes the conditional distributions\n",
    "of `alpha` for a range of values of `beta`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "betas = [10, 20, 40]\n",
    "\n",
    "for beta in betas:\n",
    "    cond = suite.Conditional(0, 1, beta)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Figure [fig.paintball3] shows the results, which we could fully describe\n",
    "as “posterior conditional marginal distributions.” Whew!\n",
    "\n",
    "If the variables were independent, the conditional distributions would\n",
    "all be the same. Since they are all different, we can tell the variables\n",
    "are dependent. For example, if we know (somehow) that `beta = 10`, the\n",
    "conditional distribution of `alpha` is fairly\n",
    "narrow. For larger values of `beta`, the distribution of\n",
    "`alpha` is wider."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Credible intervals {#credible-intervals}\n",
    "\n",
    "![Credible intervals for the coordinates of the\n",
    "opponent.](figs/paintball5.pdf)\n",
    "\n",
    "[fig.paintball5]\n",
    "\n",
    "Another way to visualize the posterior joint distribution is to compute\n",
    "credible intervals. When we looked at credible intervals in\n",
    "Section [credible], I skipped over a subtle point: for a given\n",
    "distribution, there are many intervals with the same level of\n",
    "credibility. For example, if you want a 50% credible interval, you could\n",
    "choose any set of values whose probability adds up to 50%.\n",
    "\n",
    "When the values are one-dimensional, it is most common to choose the\n",
    "**central credible interval**; for example, the central 50%\n",
    "credible interval contains all values between the 25th and 75th\n",
    "percentiles.\n",
    "\n",
    "In multiple dimensions it is less obvious what the right credible\n",
    "interval should be. The best choice might depend on context, but one\n",
    "common choice is the maximum likelihood credible interval, which\n",
    "contains the most likely values that add up to 50% (or some other\n",
    "percentage).\n",
    "\n",
    "`thinkbayes.Joint` provides a method that computes maximum\n",
    "likelihood credible intervals."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```python\n",
    "# class Joint:\n",
    "    def MaxLikeInterval(self, percentage=90):\n",
    "        interval = []\n",
    "        total = 0\n",
    "\n",
    "        t = [(prob, val) for val, prob in self.Items()]\n",
    "        t.sort(reverse=True)\n",
    "\n",
    "        for prob, val in t:\n",
    "            interval.append(val)\n",
    "            total += prob\n",
    "            if total >= percentage/100.0:\n",
    "                break\n",
    "\n",
    "        return interval\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first step is to make a list of the values in the suite, sorted in\n",
    "descending order by probability. Next we traverse the list, adding each\n",
    "value to the interval, until the total probability exceeds\n",
    "`percentage`. The result is a list of values from the suite.\n",
    "Notice that this set of values is not necessarily contiguous.\n",
    "\n",
    "To visualize the intervals, I wrote a function that “colors” each value\n",
    "according to how many intervals it appears in:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def MakeCrediblePlot(suite):\n",
    "    d = dict((pair, 0) for pair in suite.Values())\n",
    "\n",
    "    percentages = [75, 50, 25]\n",
    "    for p in percentages:\n",
    "        interval = suite.MaxLikeInterval(p)\n",
    "        for pair in interval:\n",
    "            d[pair] += 1\n",
    "\n",
    "    return d"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`d` is a dictionary that maps from each value in the suite to\n",
    "the number of intervals it appears in. The loop computes intervals for\n",
    "several percentages and modifies `d`.\n",
    "\n",
    "Figure [fig.paintball5] shows the result. The 25% credible interval is\n",
    "the darkest region near the bottom wall. For higher percentages, the\n",
    "credible interval is bigger, of course, and skewed toward the right side\n",
    "of the room."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Discussion\n",
    "\n",
    "This chapter shows that the Bayesian framework from the previous\n",
    "chapters can be extended to handle a two-dimensional parameter space.\n",
    "The only difference is that each hypothesis is represented by a tuple of\n",
    "parameters.\n",
    "\n",
    "I also presented `Joint`, which is a parent class that\n",
    "provides methods that apply to joint distributions:\n",
    "`Marginal`, `Conditional`, and\n",
    "`MakeLikeInterval`. In object-oriented terms,\n",
    "`Joint` is a mixin (see\n",
    "<http://en.wikipedia.org/wiki/Mixin>).\n",
    "\n",
    "There is a lot of new vocabulary in this chapter, so let’s review:\n",
    "\n",
    "Joint distribution:\n",
    ":   A distribution that represents all possible values in a\n",
    "    multidimensional space and their probabilities. The example in this\n",
    "    chapter is a two-dimensional space made up of the coordinates\n",
    "    `alpha` and `beta`. The joint distribution\n",
    "    represents the probability of each (`alpha`,\n",
    "    `beta`) pair.\n",
    "\n",
    "Marginal distribution:\n",
    ":   The distribution of one parameter in a joint distribution, treating\n",
    "    the other parameters as unknown. For example,\n",
    "    Figure [fig.paintball2] shows the distributions of\n",
    "    `alpha` and `beta` independently.\n",
    "\n",
    "Conditional distribution:\n",
    ":   The distribution of one parameter in a joint distribution,\n",
    "    conditioned on one or more of the other parameters.\n",
    "    Figure [fig.paintball3] several distributions for\n",
    "    `alpha`, conditioned on different values of\n",
    "    `beta`.\n",
    "\n",
    "Given the joint distribution, you can compute marginal and conditional\n",
    "distributions. With enough conditional distributions, you could\n",
    "re-create the joint distribution, at least approximately. But given the\n",
    "marginal distributions you cannot re-create the joint distribution\n",
    "because you have lost information about the dependence between\n",
    "variables.\n",
    "\n",
    "If there are $n$ possible values for each of two parameters, most\n",
    "operations on the joint distribution take time proportional to $n^2$. If\n",
    "there are $d$ parameters, run time is proportional to $n^d$, which\n",
    "quickly becomes impractical as the number of dimensions increases.\n",
    "\n",
    "If you can process a million hypotheses in a reasonable amount of time,\n",
    "you could handle two dimensions with 1000 values for each parameter, or\n",
    "three dimensions with 100 values each, or six dimensions with 10 values\n",
    "each.\n",
    "\n",
    "If you need more dimensions, or more values per dimension, there are\n",
    "optimizations you can try. I present an example in Chapter [species].\n",
    "\n",
    "You can download the code in this chapter from\n",
    "<http://thinkbayes.com/paintball.py>. For more information see\n",
    "Section [download]."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercises\n",
    "\n",
    "In our simple model, the opponent is equally likely to shoot in any\n",
    "direction. As an exercise, let’s consider improvements to this model.\n",
    "\n",
    "The analysis in this chapter suggests that a shooter is most likely to\n",
    "hit the closest wall. But in reality, if the opponent is close to a\n",
    "wall, he is unlikely to shoot at the wall because he is unlikely to see\n",
    "a target between himself and the wall.\n",
    "\n",
    "Design an improved model that takes this behavior into account. Try to\n",
    "find a model that is more realistic, but not too complicated."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
